{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a target=\"_blank\" href=\"https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%205%20-%20Scaling%2C%20Inverse%20Kinematics%2C%20and%20Inverse%20Dynamics.ipynb\">\n",
    "  <img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/>\n",
    "</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EhobipLVJg7S"
   },
   "source": [
    "# 5. Scaling, Inverse Kinematics, and Inverse Dynamics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SwX6QwGU0Jwj"
   },
   "source": [
    "## 5.1. Objectives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4p_3qlB60Np3"
   },
   "source": [
    "**Purpose**\n",
    "\n",
    "The purpose of this tutorial is to demonstrate how [OpenSim](https://opensim.stanford.edu/) [[1]](https://doi.org/10.1109/TBME.2007.901024) can be used to scale a model, and then solve the Inverse Kinematics (IK) and Inverse Dynamics (ID) problems using experimental data. \n",
    "\n",
    "To diagnose movement disorders and study human movement, biomechanists frequently ask human subjects to perform movements in a motion capture laboratory and use computational tools to analyze these movements. A common step in analyzing a movement is to compute the joint angles and joint moments of the subject during movement. OpenSim has tools for computing these quantities:\n",
    "\n",
    "*   **Inverse Kinematics (IK)** is used to compute joint angles.\n",
    "*   **Inverse Dynamics (ID)** is used to compute net joint reaction forces and net joint moments.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "urG9ZNhGHMLo"
   },
   "source": [
    "\n",
    "![Inverse Kinematics and Dynamics.png]()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6PpwC4C9HO1M"
   },
   "source": [
    "\n",
    "Inverse Kinematics computes the joint angles for a musculoskeletal model that best reproduce the motion of a subject. Inverse Dynamics then uses joint angles, angular velocities, and angular accelerations of the model, together with the experimental ground reaction forces and moments, to solve for the net reaction forces and net moments at each of the joints. The schematic above shows an overview of the inverse kinematics and inverse dynamics problems.\n",
    "\n",
    "In this tutorial, you will:\n",
    "\n",
    "*   Become familiar with OpenSim's Scale, Inverse Kinematics and Inverse Dynamics tools.\n",
    "*   Solve an Inverse Kinematics and an Inverse Dynamics problem using experimental data.\n",
    "*   Investigate the dynamic inconsistencies that arise during Inverse Dynamics.\n",
    "\n",
    "**Format**\n",
    "\n",
    "Each section of the tutorial guides you through certain tools within OpenSim and asks you to answer a few questions. The questions can be answered based on information from OpenSim and basic knowledge of the human musculoskeletal system. As you complete each section of the tutorial, feel free to explore OpenSim and the lower extremity model further on your own.\n",
    "\n",
    "**Resources**\n",
    "\n",
    "All of the files necessary to complete this tutorial are available in the following Google Drive folder: https://drive.google.com/drive/folders/1txjRDKPETE1NpOwzv0pSdIiCoc8aIZTf?usp=sharing\n",
    "\n",
    "While you complete the tutorial, you will find cells that automatically download the required files.\n",
    "\n",
    "To obtain a link to one of the files in Google Drive, right-click over it and click on **Get Link**.\n",
    "\n",
    "**Local Installation**\n",
    "\n",
    "You may also run this notebook using a locally installed OpenSim `conda` package. Visit [the OpenSim python scripting Confluence page](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53085346/Scripting+in+Python) for more information on how to install OpenSim in a local Python environment. You will also need to download the resources from Google Drive to the same directory as this notebook. Then, you may skip to Section 5.3 below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ADllcxbsjTr8"
   },
   "source": [
    "## 5.2. Set up Conda and OpenSim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mmmTuEI4HQJG"
   },
   "source": [
    "First, set up the environment by executing the following cell (See [Tutorial 1: Introduction to OpenSim in Colab](https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%201%20-%20Introduction%20to%20OpenSim.ipynb) for more information)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "0w9LYQXx5TO4"
   },
   "outputs": [],
   "source": [
    "!pip install -q condacolab\n",
    "import condacolab\n",
    "condacolab.install()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vt5loGbdfKb8"
   },
   "source": [
    "Don't worry if after executing the previous cell you get an error saying that your session has failed. This happens because `condacolab` needs to restart the session for the changes to take effect. Therefore, you have to execute the previous cell before executing any other cell.\n",
    "\n",
    "Now, execute the following cell to install the OpenSim `conda` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "vWf9ObZJ5TO5"
   },
   "outputs": [],
   "source": [
    "!conda install opensim-org::opensim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "51viyamRzzZK"
   },
   "source": [
    "## 5.3. Scaling a Musculoskeletal Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Z9iDkg-QVVhm"
   },
   "source": [
    "In this tutorial, you will usea generic musculoskeletal model with 23 degrees of freedom and actuated by 54 muscles, named _3DGaitModel2354_ [[2]](https://doi.org/10.1109/10.102791) [[3]](https://doi.org/10.1080/10255849908907988). Detailed information about the model can found on the [Gait2392 and 2354 Models](https://simtk-confluence.stanford.edu:8443/display/OpenSim/Gait+2392+and+2354+Models) OpenSim documentation page hosted on Confluence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5s5Jso_ozzZN"
   },
   "source": [
    "Subject-specific modeling involves the following steps:\n",
    "*  **Scaling** a generic musculoskeletal model to modify the body scale factors of the generic model so that it matches the anthropometry of a particular subject.\n",
    "*  **Registering** the markers placed on the model to match the locations on the subject.\n",
    "\n",
    "Scaling and registration have a significant effect on the solutions to the Inverse Kinematics and Inverse Dynamics problems, so it is crucial to perform these steps carefully."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-0nVeiJAXruu"
   },
   "source": [
    "A detailed explanation of the Scale Tool can be found on the [Scaling page](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53090000/Scaling) of the documentation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TIQp3TtGYMyQ"
   },
   "source": [
    "**Download Files**\n",
    "\n",
    "First, let's download the model file (`gait2354_simbody.osim`). The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "EAb6JtzeYObU"
   },
   "outputs": [],
   "source": [
    "!gdown \"1roNo6te-wyNk-4UPqfmV7SmUfef4Bfq3&confirm=t\" # gait2354_simbody.osim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bjlOR8fym2Pl"
   },
   "source": [
    "Load the model into OpenSim. To test that the file has been downloaded correctly, let's print its name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "oGTQE8Azh6p8"
   },
   "outputs": [],
   "source": [
    "import opensim as osim\n",
    "\n",
    "# Load the model.\n",
    "gait2354 = osim.Model('gait2354_simbody.osim')\n",
    "print(\"Name of the model:\", gait2354.getName())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-ooHO9m13zen"
   },
   "source": [
    "To scale the model, we are going to use a configuration file (`gait2354_Setup_Scale.xml`). This Scale Setup file is an XML file that contains pre-configured settings to scale and register the generic gait2354 musculoskeletal model to the dimensions of a subject for which we have experimental data. A detailed explanation of the Scale Tool can be found on the [Scaling](https://simtk-confluence.stanford.edu:8443/display/OpenSim/Scaling) page of the documentation. The following cell downloads the file. The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1ojzW5C2iHki"
   },
   "outputs": [],
   "source": [
    "!gdown \"1c0P8cN1zUmVFJxqvcpHZEDacyAoa6lQy&confirm=t\" # gait2354_Setup_Scale.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jg5Pu8k_7osp"
   },
   "source": [
    "We also need to download the MarketSet file (`gait2354_Scale_MarkerSet.xml`), which defines the position of the markers used when capturing experimental data. The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "u2MecrQliXiK"
   },
   "outputs": [],
   "source": [
    "!gdown \"1VNDVGM-MKVrYs9w3aTPCLHOfwK-s2Nn8&confirm=t\" # gait2354_Scale_MarkerSet.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wqFikw-zPau0"
   },
   "source": [
    "Finally, we need a marker file (`subject01_static.trc`) which provides information about the position of the markers during the *static trial*, which is a trial in a known pose and not moving and is a typical trial used for the scaling step. In this case, we download a file with information of the user in a static pose. The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "doLIHrRKi1KU"
   },
   "outputs": [],
   "source": [
    "!gdown \"19-ouxBQN-XUOJ-K8lRe16DUs8fqaN1Yz&confirm=t\" # subject01_static.trc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "k1bbjgNFS8h2"
   },
   "source": [
    "**Create Scale Tool**\n",
    "\n",
    "The Scale Tool helps us scale the model. We can pass a configuration file as a parameter (`gait2354_Setup_Scale.xml`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "juWpas9abaEf"
   },
   "outputs": [],
   "source": [
    "# Create a ScaleTool object using the configuration file.\n",
    "scale_tool = osim.ScaleTool('gait2354_Setup_Scale.xml')\n",
    "\n",
    "# Print some information of the config file to test everything is correct.\n",
    "print(\"Name:\", scale_tool.getName())\n",
    "print(\"Subject Mass:\", scale_tool.getSubjectMass())\n",
    "print(\"Subject Height:\", scale_tool.getSubjectHeight())\n",
    "print(\"Notes:\", scale_tool.getPropertyByName(\"notes\").toString())\n",
    "print()\n",
    "\n",
    "# Get model marker file name.\n",
    "generic_model_maker = scale_tool.getGenericModelMaker()\n",
    "print(\"Marker Set File Name:\", generic_model_maker.getMarkerSetFileName())\n",
    "print()\n",
    "\n",
    "# Get marker file name.\n",
    "marker_placer = scale_tool.getMarkerPlacer()\n",
    "print(\"Marker Placer File Name:\", marker_placer.getMarkerFileName())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6UoQfXjAJiGk"
   },
   "source": [
    "**Perform Scaling**\n",
    "\n",
    "Finally, we can scale the model by running the scale tool. You can visualize the resulting model (`subject01_simbody.osim`) in [OpenSim GUI](https://simtk.org/frs/index.php?group_id=91)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ISVbNGw-J-6_"
   },
   "outputs": [],
   "source": [
    "# Run Scale Tool.\n",
    "scale_tool.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QKEfjnyn5ldK"
   },
   "source": [
    "**Questions:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5-GGyw8X5ldQ"
   },
   "source": [
    "1. What is the mass of the generic musculoskeletal model? What was the mass of the subject? _Hint: Look for it in the `gait2354_simbody.osim` and `gait2354_Setup_Scale.xml` files._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kFTHo50B5ldQ"
   },
   "source": [
    "2. What frequency was the experimental motion data captured?\n",
    "_Hint: Look for it in the `subject01_static.trc` file._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UHMJUvVO5ldQ"
   },
   "source": [
    "3. Which body segments were scaled manually? _Hint: Look for it in the `gait2354_Setup_scale.xml` file._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ad2kW88oKPaw"
   },
   "source": [
    "## 5.4. Inverse Kinematics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sXZXzssqVmUN"
   },
   "source": [
    "Kinematics is the study of motion without considering the forces and moments that produce that motion [[4]](https://onlinelibrary.wiley.com/doi/book/10.1002/9780470549148). The purpose of Inverse Kinematics (IK) is to estimate the joint angles of a particular subject from experimental data. In this section, you will estimate a subject's joint angles during walking by performing an IK analysis using the subject scaled model and experimentally collected walking data. \n",
    "\n",
    "For each time step of recorded motion data, IK computes a set of joint angles that put the model in a configuration that \"best matches\" the experimental kinematics. OpenSim determines this \"best match\" by solving a weighted least squares optimization problem with the goal of minimizing marker error. Marker error is defined as the distance between an experimental marker and the corresponding model marker. Each marker has an associated weighting value, specifying how strongly that marker's error term should be minimized in the least squares problem. For each time step, the `InverseKinematicsTool` solves for a vector of generalized coordinates (e.g., joint angles), $q$, that minimizes the weighted sum of marker errors, which is expressed as:\n",
    "\n",
    "$\\underset{q}{\\text{min}}\\Bigg[ \\sum_{i \\in \\text{markers}}{w_i || x_i^{\\text{exp}} - x_i(q) ||^2}^{}\\Bigg] $\n",
    "\n",
    "where $q$ is the vector of generalized coordinates (e.g., joint angles), $x_i^{\\text{exp}}$ is the position of the experimental marker $i$, $x_i(q)$ is the position of the corresponding model marker $i$ (which depends on $q$), and $w_i$ is the weight associated with marker $i$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1ATo0wsJX66B"
   },
   "source": [
    "A detailed explanation of the `InverseKinematicsTool` can be found on the [Inverse Kinematics page](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53090037/Inverse+Kinematics) of the documentation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "e2kJm4w-YuMF"
   },
   "source": [
    "**Download the files**\n",
    "\n",
    "First, let's download the setup file for the `InverseKinematicsTool`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AlCIWcl6qXd2"
   },
   "source": [
    "The first file to download is the XML setup file (`subject01_Setup_IK.xml`) for the IK tool. The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "PoPq74CGj05x"
   },
   "outputs": [],
   "source": [
    "!gdown \"1z1YGOTZlTCA4jpjRFbz0FlEaen1bETou&confirm=t\" # subject01_Setup_IK.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Aw4ARIVvqhuR"
   },
   "source": [
    "The second file is the marker file (`subject01_walk1.trc`), which provides information about the position of the markers during a motion trial. The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9dJOeVIqj_aX"
   },
   "outputs": [],
   "source": [
    "!gdown \"10TVRHRNJNy766Zbejjv05n3SKh0B80qD&confirm=t\" # subject01_walk1.trc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "G6wtd-OJaF6m"
   },
   "source": [
    "**Create the `InverseKinematicsTool`**\n",
    "\n",
    "The `InverseKinematicsTool` helps us solve the Inverse Kinematics problem. We can pass a configuration file as a parameter (`subject01_Setup_IK.xml`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DaP2nPJraF6m"
   },
   "outputs": [],
   "source": [
    "# Create an IK object using the configuration file.\n",
    "inverse_kinematics_tool = osim.InverseKinematicsTool('subject01_Setup_IK.xml')\n",
    "\n",
    "# Print some information of the config file to check that everything is correct.\n",
    "print(\"Name:\", inverse_kinematics_tool.getName())\n",
    "print(\"Model File:\", inverse_kinematics_tool.get_model_file())\n",
    "print(\"Marker File:\", inverse_kinematics_tool.get_marker_file())\n",
    "print(\"Accuracy:\", inverse_kinematics_tool.get_accuracy())\n",
    "print(\"Time Range: [\", inverse_kinematics_tool.get_time_range(0), \",\", inverse_kinematics_tool.get_time_range(1), \"]\")\n",
    "print(\"Constraint Weight:\", inverse_kinematics_tool.get_constraint_weight())\n",
    "print()\n",
    "\n",
    "# Print weights information\n",
    "print(\"Weights:\")\n",
    "task_set = inverse_kinematics_tool.get_IKTaskSet()\n",
    "for i in range(task_set.getSize()):\n",
    "    task = task_set.get(i)\n",
    "    print(task.getName())\n",
    "    print(task.getWeight())\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lm09IBIac8QE"
   },
   "source": [
    "**Perform Inverse Kinematics**\n",
    "\n",
    "Finally, we run the IK tool. You can visualize the resulting motion file (`subject01_walk1_ik.mot`) by associating it to the scaled model (`subject01_simbody.osim`) in the [OpenSim GUI](https://simtk.org/frs/index.php?group_id=91)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "k6EIzdm9dEEx"
   },
   "outputs": [],
   "source": [
    "inverse_kinematics_tool.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GLp_exuTePfI"
   },
   "source": [
    "**Visualize the marker errors**\n",
    "\n",
    "The output file `subject01_ik_marker_errors.sto` provides the markers errors associated with each frame of the generated motion. You can parse this file using a `TableProcessor` and visualize the error obtained per frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aiBr4M-wxrRP"
   },
   "outputs": [],
   "source": [
    "# Use the TableProcessor to read the motion file.\n",
    "table = osim.TableProcessor(\"subject01_ik_marker_errors.sto\")\n",
    "# Process the file.\n",
    "tableErrors = table.process()\n",
    "# Print labels for each column.\n",
    "print(tableErrors.getColumnLabels())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-Ov2v2k9x5oe"
   },
   "source": [
    "We will plot the total squared error, the RMS error for the markers, and the maximum error for the markers. In the next cell you will extract this information from the motion file.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Xhnx_ulsx6CS"
   },
   "outputs": [],
   "source": [
    "# Get columns we want to plot, and the independent column (Time).\n",
    "total_squared_error = tableErrors.getDependentColumn('total_squared_error')\n",
    "marker_error_RMS = tableErrors.getDependentColumn('marker_error_RMS')\n",
    "marker_error_max = tableErrors.getDependentColumn('marker_error_max')\n",
    "x_time = tableErrors.getIndependentColumn()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WfpRt8xMjhy3"
   },
   "source": [
    "Now, we can plot the error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "w4J1bJGpzC3U"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Marker Errors from Inverse Kinematics', fontsize=16)\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(x_time, total_squared_error.to_numpy(), label='total_squared_error')\n",
    "axs.plot(x_time, marker_error_RMS.to_numpy(), label='marker_error_RMS')\n",
    "axs.plot(x_time, marker_error_max.to_numpy(), label='marker_error_max')\n",
    "axs.set_xlabel('Time (s)')\n",
    "axs.set_ylabel('Marker Error (m)')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uJ25R72HJ9YX"
   },
   "source": [
    "**Questions:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_F4VGaDkmW-I"
   },
   "source": [
    "4. Which markers have weighting values less than one? Why?\n",
    "*Hint: Think about joints that have not been modeled.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Oa7-QlzImfut"
   },
   "source": [
    "5. What is the root-mean-squared (RMS) error of all the markers in the last frame of the motion? Include units. Does this seem reasonable? Explain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hSO963Q1J7o9"
   },
   "source": [
    "6. What was the value of the maximum marker error in the last frame? Include units. Which marker had this maximum error, and why?\n",
    "*Hint: Think about the weighted least squares problem.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SzMh1MS8QazO"
   },
   "source": [
    "## 5.5. Inverse Dynamics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QL2BMH6tQiVV"
   },
   "source": [
    "Dynamics is the study of motion and the forces and moments that produce that motion [[4]](https://onlinelibrary.wiley.com/doi/book/10.1002/9780470549148). The purpose of Inverse Dynamics is to estimate the forces and moments that cause a particular motion, and its results can be used to infer how muscles are utilized in that motion. To determine these forces and moments, equations of motion for the system are solved iteratively [[5]](https://doi.org/10.1115/1.2834295). The equations of motion are derived using the kinematic description and mass properties of a musculoskeletal model. Then, using the joint angles from Inverse Kinematics and experimental ground reaction force data, the net reaction forces and net moments at each of the joints are calculated such that the dynamic equilibrium conditions and boundary conditions are satisfied [[5]](https://doi.org/10.1115/1.2834295).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cR5SdHfdYH30"
   },
   "source": [
    "A detailed explanation of the `InverseDynamicsTool` can be found on the [Inverse Dynamics page](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53090074/Inverse+Dynamics) of the OpenSim documentation hosted on confluence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wnJUuAW8RGpe"
   },
   "source": [
    "**Download the files**\n",
    "\n",
    "First, let's download the configuration file (`subject01_Setup_InverseDynamics.xml`) for the `InverseDynamicsTool`. The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "SDK4_QCkk0vo"
   },
   "outputs": [],
   "source": [
    "!gdown \"1aE_GIV0dNepL6oUGJalznsJQZDlSfcet&confirm=t\" # subject01_Setup_InverseDynamics.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XlRizvvXrXtQ"
   },
   "source": [
    "The second file is a external loads file (`subject01_walk1_grf.xml`), that contains the external loads applied to the model. The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "a3Ym-2KglEey"
   },
   "outputs": [],
   "source": [
    "!gdown \"1OZtWIpJ1UH2uViiEEeRluwfi54wsgIqG&confirm=t\" # subject01_walk1_grf.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ePpe9dGXr_II"
   },
   "source": [
    "The third file is a coordinates file (`subject01_walk1_grf.mot`), containing coordinate data. The file is stored in a Google Drive folder (See **Resources** in the Objectives section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "IoYI-ik3lUsq"
   },
   "outputs": [],
   "source": [
    "!gdown \"18L9Nj9ff-2evQkiXqD2KDZXu3DWZ30vh&confirm=t\" # subject01_walk1_grf.mot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oHeUXkItSy6D"
   },
   "source": [
    "**Create the `InverseDynamicsTool`**\n",
    "\n",
    "The `InverseDynamicsTool` helps us perform Inverse Dynamics. We can pass a configuration file as a parameter (`subject01_Setup_InverseDynamics.xml`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "r17DfremSy6E"
   },
   "outputs": [],
   "source": [
    "# Create an ID object using the configuration file.\n",
    "inverse_dynamics_tool = osim.InverseDynamicsTool('subject01_Setup_InverseDynamics.xml')\n",
    "\n",
    "# Print some information of the config file to check that everything is correct.\n",
    "print(\"Name:\", inverse_dynamics_tool.getName())\n",
    "print(\"Model File:\", inverse_dynamics_tool.getModelFileName())\n",
    "print(\"External Loads File:\", inverse_dynamics_tool.getExternalLoadsFileName())\n",
    "print(\"Coordinates File:\", inverse_dynamics_tool.getCoordinatesFileName())\n",
    "print(\"Time Range: [\", inverse_dynamics_tool.getStartTime(), \",\", inverse_dynamics_tool.getEndTime(), \"]\")\n",
    "print(\"Lowpass Cutoff Frequency:\", inverse_dynamics_tool.getLowpassCutoffFrequency())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fHgMyZc2Y6JO"
   },
   "source": [
    "**Run Inverse Dynamics**\n",
    "\n",
    "Finally, we can run Inverse Dynamics by running the `InverseDynamicsTool`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "VF3_KRo8Y6Jk"
   },
   "outputs": [],
   "source": [
    "inverse_dynamics_tool.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lUy_Vm8izR0f"
   },
   "source": [
    "**Visualize the Inverse Dynamics results**\n",
    "\n",
    "After running the Inverse Dynamics Tool, a file named `inverse_dynamics.sto` is generated. This file contains the results of the Inverse Dynamics execution, and you can parse it using a `TimeSeriesTable` to visualize the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "07K30qNozR0z"
   },
   "outputs": [],
   "source": [
    "# Load the inverse dynamics results into a TimeSeriesTable.\n",
    "tableID = osim.TimeSeriesTable(\"ResultsInverseDynamics/inverse_dynamics.sto\")\n",
    "\n",
    "# Print labels for each column.\n",
    "print(tableID.getColumnLabels())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bd3MDg8azR0z"
   },
   "source": [
    "We are going to compare the moments of the ankle angle of both legs, so let's extract these quantities from the Inverse Dynamics results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "hxG43RnlzR00"
   },
   "outputs": [],
   "source": [
    "# Get columns we want to represent, and the independent column (Time).\n",
    "ankle_angle_r_moment = tableID.getDependentColumn('ankle_angle_r_moment')\n",
    "ankle_angle_l_moment = tableID.getDependentColumn('ankle_angle_l_moment')\n",
    "x_time = tableID.getIndependentColumn()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zLngF9DkzR00"
   },
   "source": [
    "Now, we can plot the moments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ql-_ZYAizjw0"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Ankle Angle Moment Comparison')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(x_time, ankle_angle_r_moment.to_numpy(), label='ankle_angle_r_moment')\n",
    "axs.plot(x_time, ankle_angle_l_moment.to_numpy(), label='ankle_angle_l_moment')\n",
    "axs.set_xlabel('Time (s)')\n",
    "axs.set_ylabel('Moment (N-m)')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cHOJv4Mg1W_U"
   },
   "source": [
    "In solving the Inverse Dynamics problem, both kinematic data and force plate data were used, making this an over-determined problem. In other words, the problem has more equations than unknowns (i.e., degrees of freedom). Due to errors in the experimental motion data and inaccuracies in the musculoskeletal model, it turns out that Newton's second law is violated, or $\\overrightarrow{F}_{\\text{exp}} \\neq m * \\overrightarrow{a}$ [[5]](https://doi.org/10.1115/1.2834295). One method to handle this inconsistency is to compute and apply residual forces and moments to a particular body segment in the model, such that Newton's second law becomes:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZFsf78-5t2L6"
   },
   "source": [
    "$\\overrightarrow{F}_{\\text{exp}} + \\overrightarrow{F}_{\\text{residual}} = m * \\overrightarrow{a} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gFZfIUMp28Sx"
   },
   "source": [
    "Note: An analogous equation relates the ground reaction moment, $\\overrightarrow{M}_{\\text{exp}}$ to the residual moment, $\\overrightarrow{M}_{\\text{residual}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XdW9OdIA3MuC"
   },
   "source": [
    "In this musculoskeletal model, the residuals are applied to the pelvis segment. To see the residuals from the inverse dynamics solution, plot `pelvis_tx_force`, `pelvis_ty_force`, and `pelvis_tz_force` versus `time`. Using this plot, answer question 9."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_Ryg0QjYzyfU"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Get columns we want to represent, and the independent column (Time).\n",
    "pelvis_tx_force = tableID.getDependentColumn('pelvis_tx_force')\n",
    "pelvis_ty_force = tableID.getDependentColumn('pelvis_ty_force')\n",
    "pelvis_tz_force = tableID.getDependentColumn('pelvis_tz_force')\n",
    "x_time = tableID.getIndependentColumn()\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(1, 1, figsize=(7, 5))\n",
    "fig.suptitle('Pelvis Force Residuals')\n",
    "\n",
    "# Plot the knee angles on the first subplot.\n",
    "axs.plot(x_time, pelvis_tx_force.to_numpy(), label='pelvis_tx_force')\n",
    "axs.plot(x_time, pelvis_ty_force.to_numpy(), label='pelvis_ty_force')\n",
    "axs.plot(x_time, pelvis_tz_force.to_numpy(), label='pelvis_tz_force')\n",
    "axs.set_xlabel('Time')\n",
    "axs.set_ylabel('Force')\n",
    "axs.grid()\n",
    "axs.legend()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oK6v_gocKWSh"
   },
   "source": [
    "**Questions:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cC88pJIY3TLV"
   },
   "source": [
    "7. On your plot of the ankle moments, identify when heel strike, stance phase, toe off, and swing phase occur for each curve (i.e., left leg and right leg)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BLKse9bp3W7R"
   },
   "source": [
    "8. Based on your plot and the angle convention for the ankle, give an explanation of what is happening at the ankle just before toe-off.\n",
    "*Hint: It may be useful to visualize the motion in the OpenSim GUI and use the Coordinate sliders to understand the angle convention for the ankle.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nquFUQIA4E-2"
   },
   "source": [
    "9. What are the maximum magnitudes of the residual forces? Using the mass of the subject from Question 1, what fraction of body weight are the maximum residual forces?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "acVEWzHi8nRE"
   },
   "source": [
    "While applying residual forces and moments makes the model's motion dynamically consistent with the external forces, this strategy is undesirable because the residuals can be large. More advanced strategies have been developed to deal with the problem of residuals and dynamic inconsistencies, such as least-squares optimization [[5]](https://doi.org/10.1115/1.2834295), the Residual Elimination Algorithm (REA) [[6]](https://doi.org/10.1016/j.jbiomech.2005.02.010), and the Residual Reduction Algorithm (RRA) [[7]](https://bcats.stanford.edu/previous_bcats/bcats06/pdf/BCATS_2006_abstract_book.pdf#page=31). OpenSim implements a Residual Reduction Algorithm as part of its workflow for generating muscle-actuated simulations [[7]](https://bcats.stanford.edu/previous_bcats/bcats06/pdf/BCATS_2006_abstract_book.pdf#page=31). A detailed explanation of the Residual Reduction Algorithm (RRA) can be found on the Residual Reduction Algorithm page of the documentation. For additional information on these strategies, please also see [[1]](https://doi.org/10.1080/10255842.2011.627560), [[5]](https://doi.org/10.1115/1.2834295), [[6]](https://doi.org/10.1016/j.jbiomech.2005.02.010) and [[7]](https://bcats.stanford.edu/previous_bcats/bcats06/pdf/BCATS_2006_abstract_book.pdf#page=31)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "i7rhPT2id6gg"
   },
   "source": [
    "## 5.6. Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1XK9UgItWQMF"
   },
   "source": [
    "In this tutorial, you scaled a musculoskeletal model and used it to perform Inverse Kinematics and Inverse Dynamics. Then, you explored the results of each analyses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lLWgd8Noun6Y"
   },
   "source": [
    "## 5.7. Useful Links\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "m8M3dlT5xQcD"
   },
   "source": [
    "> **OpenSim Website:** https://opensim.stanford.edu/\n",
    ">\n",
    "> **OpenSim Python Scripting:** https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53085346/Scripting+in+Python\n",
    ">\n",
    "> **OpenSim API Documentation:** https://simtk.org/api_docs/opensim/api_docs/\n",
    "> \n",
    "> **OpenSim Creator Website:** https://opensimcreator.com/\n",
    "> \n",
    "> **SimTK Website:** https://simtk.org/projects/opensim\n",
    "> \n",
    "> **Biomechanics of Movement's Course:** https://www.youtube.com/channel/UCDNGy0KKNLQ-ztcL5h2Z6zA\n",
    ">\n",
    "> **Webinar on Scaling and Inverse Kinematics in OpenSim:** https://www.youtube.com/watch?v=ZG7wzvQC6eU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yIZL9H9ygQGU"
   },
   "source": [
    "## 5.8. Acknowledgments\n",
    "\n",
    "The experimental gait data were collected by Jill Higginson and Chand John in the Neuromuscular Biomechanics Lab at the University of Delaware [[8]](https://doi.org/10.1080/10255842.2011.627560). The data include marker trajectories and ground reaction forces for an adult male walking at a self-selected speed on an instrumented split-belt treadmill. Please note that the data distributed with OpenSim is from a different subject than the one described in the paper. Data collection protocols were the same for both subjects.\n",
    "\n",
    "Thanks to [OpenSimColab](https://simtk.org/projects/opencolab) project [[9]](https://doi.org/10.1080/10255842.2022.2104607) for creating the first OpenSim Conda package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2rwqw1LFq6Lz"
   },
   "source": [
    "## 5.9. References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TpbyAlHvfh8M"
   },
   "source": [
    "\n",
    "> [1].   Delp, S. L., Anderson, F. C., Arnold, A. S., Loan, P., Habib, A., John, C. T., Guendelman, E., & Thelen, D. G. (2007). **OpenSim: open-source software to create and analyze dynamic simulations of movement.** *IEEE Transactions on Bio-Medical Engineering*, 54(11), 1940–1950. https://doi.org/10.1109/TBME.2007.901024\n",
    "> \n",
    "> [2].   Delp, S. L., Loan, J. P., Hoy, M. G., Zajac, F. E., Topp, E. L., & Rosen, J. M. (1990). **An interactive graphics-based model of the lower extremity to study orthopaedic surgical procedures.** *IEEE Transactions on Bio-Medical Engineering*, 37(8), 757–767. https://doi.org/10.1109/10.102791\n",
    ">\n",
    "> [3]. Anderson, F. C., & Pandy, M. G. (1999). **A dynamic optimization solution for vertical jumping in three dimensions.** *Computer Methods in Biomechanics and Biomedical Engineering*, 2(3), 201–231. https://doi.org/10.1080/10255849908907988\n",
    ">\n",
    "> [4]. Winter, D. A. (1990). **The biomechanics and motor control of human movement** (2a ed.). *John Wiley & Sons*. https://onlinelibrary.wiley.com/doi/book/10.1002/9780470549148\n",
    ">\n",
    "> [5]. Kuo, A. D. (1998). **A least-squares estimation approach to improving the precision of inverse dynamics computations.** *Journal of Biomechanical Engineering*, 120(1), 148–159. https://doi.org/10.1115/1.2834295\n",
    ">\n",
    "> [6]. Thelen, D. G., & Anderson, F. C. (2006). **Using computed muscle control to generate forward dynamic simulations of human walking from experimental data.** *Journal of Biomechanics*, 39(6), 1107–1115. https://doi.org/10.1016/j.jbiomech.2005.02.010\n",
    ">\n",
    "> [7]. John, C.T., Anderson, F.C., Guendelman, E., Arnold, A.S., Delp, S.L. (21st October 2006). **An algorithm for generating muscle-actuated simulations of long-duration movements.** *Biomedical Computation at Stanford (BCATS) Symposium*, Stanford University, Poster Presentation. https://bcats.stanford.edu/previous_bcats/bcats06/pdf/BCATS_2006_abstract_book.pdf#page=31\n",
    ">\n",
    "> [8]. John, C. T., Anderson, F. C., Higginson, J. S., & Delp, S. L. (2013). **Stabilisation of walking by intrinsic muscle properties revealed in a three-dimensional muscle-driven simulation.** *Computer Methods in Biomechanics and Biomedical Engineering*, 16(4), 451–462. https://doi.org/10.1080/10255842.2011.627560\n",
    ">\n",
    "> [9] Mokhtarzadeh, H., Jiang, F., Zhao, S., & Malekipour, F. (2022). **OpenColab project: OpenSim in Google colaboratory to explore biomechanics on the web.** *Computer Methods in Biomechanics and Biomedical Engineering*, 1–9. https://doi.org/10.1080/10255842.2022.2104607"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "opensim",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
